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ABSTRACT 


Methods for obtaining approximate solutions for the fundamental eigenvalue 
of the Laplace-Beltrami operator (also referred to as the membrane eigenvalue 
problem for the vibration equation) on the unit spherical surface are developed. 
Two specific types of spherical surface domains are considered: (1) the interior 
of a spherical triangle, i.e., the region bounded by arcs of three great circles, 
and (2) the exterior of a great circle arc extending for less than tt radians on 
the sphere (a spherical surface with a slit). In both cases, zero boundary con- 
ditions are imposed. In order to solve the resulting second-order elliptic par- 
tial differential equations in two independent variables, a finite difference 
approximation is derived. The symmetric (generally five-point) finite difference 
equations that develop are written in matrix form and then solved by the iterative 
method of point successive overrelaxation. Upon convergence of this iterative 
method, the fundamental eigenvalue is approximated by iteration utilizing the 
power method as applied to the finite Rayleigh quotient. The implementation of 
these numerical techniques is described in detail. Including the presentation of 
separate algorithms for calculating the fundamental eigenvalue for the two dis- 
tinct domains considered. Although analytical solutions to this eigenvalue prob- 
lem are not available in the general case, the problem is solved analytically for 
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exact values of the fundamental eigenvalue for certain special cases. These 
exact solutions are useful in providing checks on the numerical results pre- 
sented based upon digital computer calculations. Several tables of numerical 
applications for various cases of interest are displayed and discussed. 
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SECTION I 


INTRODUCTION 

Let G be a polyhedron in Euclidean three-dimensional space whose boundary 
is denoted by BG. Consider the problem: 

Av = f in G, 

( 1 ) 

V = 0 on 3G, 


where A denotes the Laplacian and v and f are defined in G. It is known 
(Reference 1) that the derivatives of v may become singular near a vertex of 
G, and the severity of the singularity is determined by the fundamental eigen- 
value of an associated eigenvalue problem for the Laplace-Beltrami operator 
on the unit sphere. Methods for obtaining approximate solutions of this asso- 
ciated eigenvalue problem will be considered in this paper. 

If spherical co-ordinates are introduced such that 

X = r s in 0 cos 6 ^ 

y - X sin 0 sin (2) 

z - r cos 0* 

then the Laplacian may be written (Reference 2, p. 225) 


where 


Au - r“2(r2y ) + r ^Au, 

V r ' r 


(3) 


Au r CSC 0[(u^ CSC 0)^ + (u^ sin 0)^ , ^ 

and where the subscript notation is used to indicate partial differentiation. If 
the origin of the co-ordinate system is placed at a vertex of G, then the singu- 
larity in the solution of problem (1) at the origin is related to the eigenvalue 
problem 


Au -I- Au - 0 in D, 
u = 0 on 3D, 


( 5 ) 


1 
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where D represents the area on the surface of a small sphere (centered at the 
origin) and bounded by the polyhedron G. In this paper, the case in which the 
region D is a spherical triangle T, i,e., the region bounded by three great 
circles, will be considered. Also to be discussed is the case in which D is a 
slit domain, consisting of the exterior of an arc of a great circle on the sphere. 

Without loss of generality, it may be assumed that the sphere is of unit 
radius, so that the spherical co-ordinates satisfy 

r = 1, 

0 ^ ^ 77 , (6) 

0 ^ 9 < 2t7. 

The boundary 3D will consist of arcs of three great circles on the unit sphere. 
The spherical co-ordinate system is defined so that the origin of co-ordinates 
is at the center of the sphere, the z-axis intersects the unit sphere at a vertex 
of T (i.e., at an intersection point of two of the great circles specifying T), and 
the X, z -plane contains a side of T which is less than n radians in length. In 
this manner, the relevant arcs of the three great circles specifying T are given 


6-0, 

(7a) 

0 = 0, 

(7b) 

z = ax + by, 

(7c) 


where 0 , a, and b are all constants (see Figure 1). Equation (7c) may be trans- 
formed to spherical co-ordinates by use of equations (2) as 

cot 4> = a cos ^ + b sin 6. (8) 

Equations (7a), (7b), and (8) may be regarded as defining two domains on the 
sphere, the region R(0 , a, b) given by the inequalities 




led by arcs of three great circles 
^hich are specified by the three 
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0 < d < 

cot 0 > a cos (9 + b sin 6 , 


( 9 ) 


and the region S{0 , a, b) complementary to R<© , a, b). Both the interior R and 
the exterior S of the spherical triangle T are thus specified by the three param- 
eters 0 , a, and b, which satisfy the inequalities 


0 ^ 0 < 2tt, 


- 00 < a, b < 00 . 

The eigenvalue problem (5) for the unit sphere may be written, by virtue of 
equation (4), as 


( 10 ) 


(u^ CSC cf>)^ + (u^ sin 0 )^ + \u sin 0 = 0 in D, 
u = 0 on 3D, 

where D is one of the regions R(0 , a, b) or S(0, a, b). Equation (11a) may be 
expanded as 

CSC 0 + (u^ sin 0)^ + \u sin 0 = 0. 

It is well known (Reference 2, p. 298) that the eigenvalue problem (5) has a 
denumerably infinite sequence of positive eigenvalues which may be ordered so 
that 


0 < ^ ^3 i 

as well as a corresponding sequence of linearly independent eigenfunctions u^, 
Uj, . . . The smallest positive eigenvalue for this problem is known as the 
fundamental eigenvalue. In this paper, the fundamental eigenvalue of equation 
(12) will be determined for the following two cases: (1) when D is the interior 
R of a spherical triangle, and (2) when D is the exterior S of a spherical 
"triangle" which has degenerated to a line in the special situation: 0 = 0. In 
both cases, the boundary condition (11b) applies. In order to solve the 
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second-order elliptic partial differential equation (12) in two independent vari- 
ables 8 , 4>, a. finite difference approximation is derived. The resulting finite 
difference equations are written in matrix form and then solved by the iterative 
method of point successive overrelaxation. Upon convergence of this iterative 
method, the fundamental eigenvalue is found by iteration utilizing the power 
method as applied to the finite Rayleigh quotient. Numerical results for a num- 
ber of cases are presented and then discussed. 



SECTION II 


FINITE DIFFERENCE APPROXIMATION 
In the following, the domain D on the spherical surface is taken to be the 
interior R(© , a, b) of a spherical triangle, with 0 < 0 < 2tt. Then the two- 
dimensional spherical surface may be mapped directly onto the ^,0 -plane in the 
manner that cartographers would refer to as Mercator's projection in the case of 
the Earth's surface. The spherical wedge-shaped surface bounded by the merid- 
ional arcs i9 = 0, 0 is thus mapped into the plane rectangle boimded by the two 
pairs of parallel lines, 0 = 0,0 and 0 = 0,7r . The arc of the great circle form- 
ing the third side of the spherical triangle, specified by equation (8), is mapped 
into a continuous curve extending from the ordinate axis 0 = 0 to the parallel 
side of the rectangular region, 0 =tt , and lying wholly within the rectangular 
closed region (see Figure 2). 

In problems involving elliptic partial differential equations for which an 
analytic solution is not known, such as equation (11), the method of finite differ- 
ences is most commonly employed to determine numerical results. This finite- 
difference technique involves first establishing a network of grid points through- 
out the rectangular region of interest occupied by the independent variables d 
and0. Suppose, for the moment, that positive constant grid spacings of h and 
h^ are chosen in the 0- and 0 -directions, respectively. Then the rectangular 
network consists of the grid points 

( 0 ,. 0.) = (ih^, jh^), i = 0, 1. 2 N^; 

i = 0. 1, 2. N^. 

Here the positive integers N^ , N^ represent the number of grid intervals (or 
grid points less one) in the 0- and 0-directions, respectively, so that 
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<t> 



Figure 2. Network of grid points d .^ , <p. superposed over the spherical 
wedge-shaped surface containing the spherical triangle mapped onto a 
plane rectangular region 
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In this manner, the continuous problem (11) is discretized by the replacement of 
the connected open domain R of the independent variables 6 , 4 > with a finite set 
of grid points ft , (see Figure 2). If the exact solution to the partial differ- 
ential equation (12) is denoted u = u(0,<^), then let its approximation, to be de- 
termined at each grid point by the method of finite differences, be 
U = U(ft- , !^> ) = U. . . The partial derivatives involved in equation (12) will then 
be approximated by suitable finite-difference expressions involving h^ , and 
Uj j . This procedure leads to a finite system of simultaneous algebraic equa- 
tions in the Ui , j , whose values may then be determined. 

In order to approximate the second-order partial derivative in the first 
term of equation (12), use is made of the centered second difference quotient, 


which has an error of order h|. Formula (15) arises by performing a Taylor's 
series expansion about the central value U. . (Reference 3, pp. 430-431). In 
considering partial derivatives in the 4 > -direction, the nature of the singularity 
at the north pole, where two boundaries of the spherical triangle meet at a point 
with a discontinuous tangent (i.e., a corner), must be taken into account. For 
this reason, variable grid spacing in the (^-direction is introduced. Define 


^.j >0, j = l. 2, , (16) 

as the variable grid spacing parameter in the 0 -direction. Then equations (13) 
and (14b) are seen to hold only for the case of constant grid spacing in the 
0-direction. Now further define the midway (or averaged) grid points, 




'^3-1/2 +‘^3-1^ 


h " 1 . 2 , 


• - 1 - 


(17) 



9 


With these further definitions, the arrangement of grid points is as shown in Fig- 
ure 3. The first-order partial derivative appearing in the second term of equa- 
tion (12) may now be approximated by a centered midway first difference quotient. 


sin 0)^ 


"T j + 1 ^ j ^ 


(18) 


In equation (18), the error is of order h^ . , and in order to approximate the 
terms involving u^, a centered first difference quotient is employed: 


u. - u. . 

/u 'I = 

<#’^i,j + l/2 h 

j + 1 


j-1/2 


u. . - u. . , 

l.J l.J-1 

• 

<p, 1 


(19) 


Combining equations (18) and (19) results in the approximation. 


2 sin0.^, U. - U. . 
(u^ sin 0)j, = •* ‘ 


^ h, + h, . 

<p. j + 1 <p, : 


j., 

</>. j + i 


2 sin0 U. . - U. . , 

- 1/ 2 . i.j i.J-l 


^<f>, j + i ^0. j 


h , . 

1 


( 20 ) 


The finite difference approximation to equation (12) may now be written, by use 
of equations (15) and (20), as 


- 2U; 


h| \~..j + i / 




Vi + i + 


<^. j 


+ /VU- J s in 0. =0. 


( 21 ) 




• GRID POINTS 



• MIDWAY (AVERAGED) 
GRID POINTS 


Figure 3. Arrangement of grid points with grid 
spacing parameters indicated 
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Multiplying equation (21) by the factor (h^ ^ ^ j ) re-arranging terms 

yields a symmetric five-point difference equation of the form: 


where 


a. U. . - b.U. • - c.U. - b.U . • - c. U. . 

j i.j j i+i.j 1 i.j+i j i-i.j j-i i.j-i 

= \e.U. 


CSC 4>. 

b. = 1 (h. . . , + h , 


( 22 ) 


(23a) 


and 


c, 

J 


^0. j+1 


5. - (h , . . , + h , , ) s i n 0. , 

J ' ip, J + 1 <p. J ' ] 


(23b) 

(23c) 


a. = 2b. + c. + c..^. 


(23d) 


In terms of grid points in the 0 -direction, by virtue of equations (16) and (17), 


<6 - 

b ^ _}JJ izl > 0, (24a) 

^ h| s i n 0. 


and 


c, 

J 



ini(0. +0.,^) 

ci. - <i. 

^ j + 1 


> 0 , 


(24b) 


- 4>- ,) sin 0. > 0. 

) '^j + l 

It is to be noted that equations (24a) and (24c) hold for j = 1,2 “ 1* while 

equation (24b) holds for the less restrictive range, j = 0, 1, 2 - 1. In 

this way, equation (23d) also is valid for j = 1, 2 - 1. The finite dif- 

ference equation (22) is an approximation to the continuous partial differential 
equation (12), which is valid for the interior R of the spherical triangle. Thus, 
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equation (22) holds for all rectangular interior grid points, i = 1, 2, . . . , - 1; 

j = 1, 2, . . . , - 1- (Grid points exterior to the spherical triangle southern 

boundary arc will be taken into account in the following discussion.) 

The boundary conditions (lib) must be imposed on the discretized version 
(22) of the eigenvalue equation. With reference to Figure 2, it is seen that 

for j = 0, 1, 2 , (25) 

along the meridional boundaries (7a) and (7b) , respectively, and U. . = 0 for 
those grid points (i,j) such that 

cot 0. a cos + b sin (9^. (26) 

In the relation (26), equality occurs for grid points that coincide with the arc of 
the great circle forming the southern boundary (8) of the spherical triangle, while 
inequality occurs for grid points in the exterior of the triangle. Furthermore, 

for i= 0, 1, 2, ... , (27) 

at the north and south poles of the sphere, respectively, since the former lies on 
the boundary of the spherical triangle and the latter in the exterior of the triangle. 
In the limiting cases of a, b (either or both), the triangle boundary, accord- 

ing to equation (8), includes the point 0 = v, the south pole, and the boundary con- 
dition (27) is still valid. 

Variable grid spacing was adopted in equation (16) for the 0-direction in 
order to allow for the effect of the singularity at 0= 0, the north pole. In accord 
with this choice, grid points in the 0 -direction are specified by 
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2}Yn 


— i) — for j = 0, 1, 

2 2 


forj = ~N^+ 1. _N^+ 2,..., N^, 


where is assumed to be an even integer and 7 is a positive constant param- 
eter which defines the <;i>-grid spacing. If 7= 1 , then formula ( 28 ) reduces to a 
uniform grid spacing in the ^6-direction, as specified in the network formulas 
( 13 ) and ( 14 b) and for which the integer need not be even. If 7 > 1 , the 
density of grid points increases toward the poles, 4 > r . - 0 and 0 =77 , and de- 

creases toward the equator, ^6^ /, " Conversely, if 7 < 1, then the density 

4 > 

of grid points increases toward the center and decreases toward the poles. For 
all 7, formula ( 28 ) specifies grid points in the 0 -direction which are symmetri- 


cally placed about the center point or equator. 



SECTION m 


ITERATIVE SOLUTION BY SUCCESSIVE OVERRELAXATION 
In order to apply the iterative method of solution by point successive over- 
relaxation to the finite difference equation (22), it is advantageous to adopt a 
matrix formulation of the problem. Let U represent the ordered set of unknown 
values U. . written as a vector of (N^ - 1)(N^ - 1) components defined on the 
network of rectangular interior grid points (i9. N^ 1 i 

j = 1, 2, 3, . . . , N^ - 1. Then the finite difference equation (22), with the bound- 
ary conditions (25), (26), and (27) included, may be written in matrix form as 

AU = \EU. (29) 

where A and E are square matrices of order (N^ - 1)(N^ - 1). 

An example of this matrix formulation for the special case of N^ = 4 , N^ = 6 
appears in Figure 4. The vector U has a component associated with each rec- 
tangular interior grid point, and these components are numbered by a so-called 
natural ordering of the grid points (Reference 4, p. 187; Reference 5, p. 454) in 
which successive horizontal grid lines (see Figure 2) are scanned in order of 
increasing i index with the j index increasing on each new line, rectangular 
boundary grid points excluded. In this way, the r*^ component of U is associ- 
ated with the grid point , 4 >-^ ) such that r = i + (j - 1)(N^ - 1). 

In general, the matrix A will be sparse, with the non-zero elements a. 
along the principal diagonal, -b . along the two adjacent diagonals (with zeroes 
interspersed along these adjacent diagonals), and -c. along two other diagonals 
symmetrically removed from the principal diagonal by (N^ - 1) elements, the 
number of rectangular interior grid points along a horizontal grid line. Since 
Sj , bj , and Cj are all positive, by equations (23d) and (24), A consists of posi- 
tive diagonal entries and non-positive off-diagonal entries. Furthermore, it is 
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15 



u 

E 

Figure 4. The square matrix A of order (N^ - 1)(N^ - 1) = 15 with 
tridiagonal and diagonal submatrices indicated, the vector U of 15 
components (the superscript "T" indicates transpose), and the 
square diagonal matrix E of order 15 for the special case = 4, 


= [U,,U,.U„U.,U„U„U„U,,U„U,.U,.U,,U,.U,.U, 


di ag [e e e e e e e e e e e e e^] 
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seen that A is a symmetric matrix. There is a partitioning of A, as indicated 
in Figure 4, which results in the following block tridiagonal form: 



The partitioning ofAinthis manner results from placing all rectangular interior grid 
points along a particular horizontal grid line into a block. In equation (30) , each 

square diagonal submatrix A.,j=l,2 N^-lis itself tridiagonal and of 

order (N^ - 1). The off-diagonal submatrices Cj , j = 1, 2, . . . , - 2 are 

likewise square and of order (Ng, - 1) but diagonal in form. These character- 
istics follow from the fact that equation (22), which determines the form of A, 
is a five-point symmetric difference approximation. The non-zero entries of 
the submatrix A. consist of the coupling coefficients of a grid point with its 
immediate grid neighboring points on the same horizontal line, while the sub- 
matrix Cj consists of diagonal entries representing the coupling coefficients 
of a grid point with its immediate grid neighboring points on the horizontal lines 
vertically above and below. 

The matrix E in equation (29) is square, of the same order as A, and diag- 
onal with positive entries, by virtue of equation (24c). Finally, note that the 
entries of A and E depend only on the j-index and are independent of i. 

Now, define a square matrix V of order (N^ - 1)(N^ - 1) by V = 
where E ^ is the diagonal matrix consisting of entries which are equal, re- 
spectively, to the square roots of the elements of E. With this definition, equa- 
tion (29) may be transformed to 
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XV= E'l/2Arl/2Vr:XV, (31) 

through pre-multiplication by the diagonal matrix E , the inverse of 
formed by taking reciprocals of the respective diagonal entries of E The 
matrix A = , transformed from A by the pre- and post- 

multiplication of a diagonal matrix, retains the symmetry property of A. 

In order to determine the fundamental (or smallest) eigenvalue of the matrix 
X, an iterative procedure known as the power method (Reference 5, pp. 147 ff; 
Reference 6, pp. 355-356) will be utilized. As applied to problem (31), the 
algorithm for the power method involves the sequence of vectors defined by 

m=l, 2, ..... (32) 

where the sequence of scalars is given by 


;^(m) ^ (V<">), yC”)) ^ 


m = 1, 2, . . . . 


(33) 


The limit of the vectors as m-* is the eigenvector associated with the 

fundamental eigenvalue of The fundamental eigenvalue, in turn, is approximated 

by the finite Rayleigh quotient (33), which involves a ratio of inner products, in 
terms of the original vector U of equation (29), it is seen that ~ E 
and A”^ - E^^^ A”i E so that equation (32) becomes 


U(">+i) = M"'>A‘»EU(“\ m= 1, 2,.... 

Furthermore, the Rayleigh quotient (33) in terms of the original vector U is 
equivalent to 


(34) 


^(m) ^ (uC"), EUC”>) 

(A'iEU('">, EU^™)) 


m = 1, 2, . ... 


(35) 


As an initial estimate for the iterations (34), is defined so that components 

associated with grid points interior to the spherical triangle T (a subset of the 
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rectangular interior grid points enumerated by the subscripts i=l,2,3 

- 1; j = 1, 2, 3, .... - 1) are assigned the value unity and components 

associated with grid points coincident with the spherical triangle boundary or 
exterior to T are assigned zero values. That is, u/ ^ those rectangular 

interior grid points satisfying relation (26), and 11^? =1 for those rectangular 
Interior grid points which satisfy the reverse inequality; 


cot > a cos 9 .^ + b sin 9 ^. 

In order to implement the power method as expressed by equations (34) and 
(35), it is necessary to solve equations of the general form AU = F for the vector 
U, where F represents the vector iterates . For this purpose, the iterative 

method of point successive overrelaxation is applied to the system of linear equa- 
tions (22). This method (Reference 4, pp. 58-59, for example) produces iterates 
whose components are given by 




^(k+i) _ CO 

a. 

j 

+ k=l, 2 ,.... <36) 

In equation (36), w is the scalar relaxation factor (evaluated by a method to be 
described), and the superscripts k and (k+1) indicate an Iteration process dis- 
tinct from the m-iterates specified in the power method described above. In fact, 
the vectors U and F that appear in equation (36) depend upon the iteration index 
m, although this has not been indicated explicitly, in order to preserve legibility. 
Note that, according to equations (23d) and (24), a^ ¥ 0 , so that singularity prob- 
lems do not arise. The k-iterations arising from use of equation (36) will be 
called the inner iterations, so as to distinguish them from the m-iterations 
arising from use of equations (34) and (35), which will be called the outer itera- 
tions. The vectors F are defined by F^"*^ = and are k-independent 



19 


(i.e., constant for the inner iterations). Since E = diag (e^ ), the following compo- 
nent relationship holds: 


Fi_j=e.U.j, fori = 1,2 - 1; j - 1, 2,- - • , - 1. (37) 

The range of indexing for the i, j subscripts shown in equation (37) also applies 
to the inner iterations (36). During the inner iterations, the convergence 
parameters, 


ef*) = max 

1 ^ i <N^ - 1 




U 


,(k) 


k = 1. 2,... 


(38a) 


and 




k = 2, 


3,.... 


(38b) 


are calculated for use in the evaluation of the relaxation factor co . Before dis- 
cussing this evaluation procedure, it is worthwhile to note that the theoretical 
demonstration of convergence of the successive overrelaxation iterates (36) for 
any initial vector chosen follows from the Ostrowski-Keich theorem 

(Reference 7, p. 123) for a relaxation factor in the range 0< < 2. This 

theorem depends upon the positive definite property of the matrix A, which in 
turn follows (Reference 4, p. 23) from the fact that A may be shown to be a 
symmetric irreducibly diagonally dominant matrix with positive diagonal entries, 
by virtue of equation (23d) and certain graph theoretic arguments (Reference 4, 
p. 20). 

Although the method of successive overrelaxation will converge for all re- 
laxation factors such that 0 < a> < 2, the most rapid convergence occurs for 
an optimal value denoted where 1 < < 2. A theoretical expression 

for exists (Reference 4, p. 110) since the Jacobi matrix associated with 

the matrix A is cyclic of index 2. This expression depends on the spectral 
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radius of the associated Jacobi matrix, a quantity that is not, unfortunately, known 
a priori. In order to overcome this difficulty in determining , described 
{Reference 6, p. 257) as "perhaps the most important problem" in the practical 
use of successive overrelaxation, a numerical technique (Reference 6, pp. 369-370) 
utilizing the point Gauss-Seidel iterative method is applied. If co is set equal to 
unity in equation (36), then the second term on the right side of the equation van- 
ishes and the successive overrelaxation method reduces to the point Gauss- 
Seidel method. In this case, the ratio of maximum component norms for succes- 
sive inner iterates of the vector U defined in equations (38) will converge to a 
value, 


lim = r i 1, 
k-oo 

equal to the square of the desired spectral radius of the associated Jacobi matrix 


mentioned above. Then can be computed by the formula 

2 


0) 


op t 


(40) 


1 + /I - r 

A method for defining a vector relaxation factor a> =co. , dependent upon the grid 
spacing parameters and 0j in the <^-direction and not requiring Gauss-Seidel 
iterations, for use in the successive overrelaxation method (36) was also 
attempted, but without convergence success. 

In the implementation of the numerical technique described above, four addi- 
tional convergence parameters are selected initially before the finite difference 
approximation is applied. Two of these parameters, and e .are 

specified small positive numbers such that 


® ^ ^INNER’ ^OUTER ^ 


(41) 


and which are used as criteria for evaluating the convergence of the inner and 
outer iteration schemes, respectively. Two other parameters, and n\„ax * 
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are specified positive integers, the first of which is used to effect the convergence 
of r in equation (39) and also as an upper bound on the number of inner iterations 
permitted and the second of which is used simply as an upper bound on the num- 
ber of outer iterations permitted. These four pre-selected parameters are utilized 
in the numerical implementation as follows, in a modification of a method sug- 
gested (Reference 6, pp, 375-376) for the solution of the finite eigenvalue problem 
on a digital computer. First, the point Gauss-Seidel iterative method (36), in 
which = 1, is applied using the initial vector estimate , as previously de- 
scribed. After each such inner iteration is completed, the norm e is com- 
puted by equation (38a). If for some value of k in the range 1 ^ k £ k^^^ (say 
k= K), it is found that , then the inner iterations are regarded as 

having converged and r is taken to be the value r^*^^ and is computed by 

formula (40). Henceforth, all subsequent inner iterations use the point successive 
overrelaxation method (36), in which co ~ (r(*^> ), and an initial vector esti- 
mate is computed as described below. If, however, eC' > > for all 

k = 1, 2, . . . , k,„ax , then the inner iterations do not converge and r is set equal 
to r and is computed using this value for r. Henceforth, all subse- 

(k ) 

quent inner iterations use successive overrelaxation in which o) (r ). 

In the latter instance, the non-convergence of the inner iterations is disregarded 
during the first (m « 1) outer iteration. However, if, in any subsequent (m = 2, 

3, . . . ) outer iteration, the inner iterations fail to converge (i.e., ^inner 

for all k = 1, 2, . . . , k^g,^), then the eigenvalue problem is regarded as non- 
convergent. In such case of problem non-convergence, the problem may be 
attempted again by properly adjusting the pre-selected convergence parameters, 
e.g., by increasing the value of k^^^ or that of or both. 

When convergence of the inner iterations is achieved (or, alternatively, after 
k inner iterations during the first outer iteration), then the Raylei^ quotient 

max 
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(35) may be used to approximate the fundamental eigenvalue for the m*^^ 
outer iteration. The vector inner products on the right side of equation (35) may 
be written in component form as = P/Q, where 






and 


i=l \j = l 


Nfl-1 /Nx-1 


Q- (A’‘EU<‘”>, EU^"-)) = ^ f ^ F. .\}. 


(42a) 


(42b) 


In the double summation expression (42a), the vector EU is given by the compo- 
nents F. . f while the vector U (prior to initiation of the most recent set of inner 
iterations) must be retrieved by the ratio j j of components. In the double 
summation expression (42b), the vector A * EU (following completion of the most 
recent set of inner iterations) is given by the components Uj j upon convergence 
of the inner iterations. Now, for the second and all subsequent outer iterations, 
convergence of the outer iterations and of the approximation to the fundamental 
eigenvalue occurs when 

^<">-^<"-‘>1 < w (43, 

In this event, A.^’"^ is the calculated fundamental eigenvalue which approximately 
satisfies the eigenvalue problem (11). If, however. 


( 44 ) 

for all m = 2, 3, . . . , m^^^ , then the fxmdamental eigenvalue does not converge 
within the allotted number m of outer iterations. This situation may be 
alleviated by properly adjusting the pre-selected convergence parameters, e.g., 
by increasing the value of m^^^ or that of ^quter both. Finally, if inequality 
(44) is satisfied for some outer iteration m, where 2 s m < m„a,t , then the outer 
iteration procedure has not converged, but the vector U is re-initialized for the 
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inner iterations utilizing point successive overrelaxation (36) by the power method 
(34) , which may be written in component form as 


^ for i = 1, 2,-... - 1; j = 1. 2,..., N^- 1. (^5) 

In equation (45), the vector components on the right side are those that result 
following completion of the inner iterations for outer iteration m. The vector 
components on the left are those to be utilized as the initial vector estimate for 
the first (k = 1) inner iteration of outer iteration number (m + 1) as well as the 
components to be used in re-evaluating the vector F by equation (37). 

This completes the discussion of the iterative method of solution of the finite 
difference equations by point successive overrelaxation and the subsequent iter- 
ative determination of the fundamental eigenvalue by the power method utilizing 
the Rayleigh quotient. Application of these techniques is clarified by the algo- 
rithmic presentation within Appendix A of the methods discussed in Sections II 


and m. 



SECTION IV 


ANALYTIC SOLUTIONS FOR SPECIAL CASES 


In this section, the elliptic partial differential equation (12) will be solved 
analytically, or exactly, for certain special cases, although, as has been re- 
marked previously, an analytic solution is not known for the general case. The 
exact solutions to be found will be useful in providing checks upon numerical 
calculations based upon the finite difference approximation. 

In order to solve equation (12) analytically, the separation of variables 
technique will be employed (Reference 2, pp. 510-512), in which 


u(^, 4>) ^ X(6) Y(4>)- 

Substitution of the separated form (46) into equation (12) yields, after 
re-arrangement , 


1 d2X((9) _ sin (p 

X(^) d^2 " *yW 


d^Y(<ifc) 


sin (p ^ 


dY(0) 


d(p 


cos 4>]+ X s i n2 (p. 


(47) 


Since the left side of equation (47) is a function of d only, and the right side of 
the equation is a function of 4> only, both sides may be set equal to a separation 
parameter k , a constant independent of both 9 and <p. In this case, the equation 
for 9 , 


+ k^X(i9) = 0, 

69^ 

is readily seen to have the solution 

X(9) = Asink(9 + B cos kd, (48) 

where A and B are constants. The boundary condition (lib) requires that 
u(O,0) = u(0 ,<p) = 0, so that the solution (48) becomes 


XC^) = A sin 



(49) 
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where k=s7r/0 is chosen so that X(0) vanishes only at the endpoints of the in- 
terval 0 £, & £. 0 . The equation for 4>i from equation (47) , 

s in^ ^ ^ +sinc^cos0 + (X. sin^ 4> - Y(0) 0, 

d02 <30 

may be transformed by the substitutions cos 4>= Y(cos 0) *= v(z) to the follow- 
ing form: 

(1 - 22)2 £2lf) _ 2z(l - z2) + [\(1 - z2) - k2] V(z) r Q. (50) 

dz2 

The solutions to equation (50) are known (Reference 2, pp. 327, 511) to be the 
associated Legendre functions of order k, 

= Pn.k(z) = Cl - — P„(z). k = 0, 1, 2,..., n, (51) 

dz*< 

where n is defined by X = n(n+ 1). The Legendre functions of order zero, 

P (z) = P ^ (z), are polynomials in z of degree n. 

Consider the special case k = n. For this choice, the derivative factor in 
equation (51) reduces to a constant and v(z) becomes a multiple of (sin0)’‘. In 
this case, the solution (46) is 

u(0. <jb) = Asin^^j(sin<^)’^, (52) 

where A is an arbitrary constant, and the corresponding eigenvalue is 



Note that for A^O, the eigenfunction (52) vanishes only for (9 = 0,0 and 0 = 0, v 
in the domain 0 < 9 < 0 < 27 t; 0 <</> < tt* Thus, equation (52) represents a solu- 
tion for the case in which the spherical triangle has degenerated to a two-sided 
spherical wedge-shaped domain bounded by meridional arcs extending from 

the north to the south pole» In fact, it may be readily demonstrated that the 



eigenfunction (52) satisfies equation (12) by direct substitution, and, as noted, 
it further satisfies the boundary conditions (11b). Hence, the eigenfunction (52) 
is an exact anal 3 i;ic solution to problem (11), where 



for any 0 in the range 0 < 0 < 2tt. The original restriction of the solution (51) 
that k = n = 7 t/ 0 be a positive integer may thus be removed in this special case. 
In terms of the constant parameters a,b of inequality (9), the wedge-shaped 
domain is specified by the limits a = 0,b when 0 < tt ; however, when 0 > 77, 
the wedge does not appear as a limiting case of the spherical triangle. 

Now consider the special case k = n - 1. For this choice, the derivative 
factor in the associated Legendre function (51) reduces to a term linear in z. 
This term linear in z is, in fact, merely a constant multiple of z, since the 
Legendre functions P„ (z) are polynomials in z containing terms only of even 
or odd powers of z, depending as n is even or odd, respectively. It is seen 
then that v(z) is a multiple of (sin 0)*' cos 4 > i° this case, and the solution (46) 
is 

u(^, 0) = Asin^^j(sin0)^/® cos<;6, (53) 

where A is an arbitrary constant, and the corresponding eigenvalue is 

X = n(n + 1) - (k + 1) (k + 2) = (.^ + ij + 2 

For A ^ 0, the eigenfunction (53) has the same zeroes as that of the previous 
eigenfunction (52), viz., (9 = 0, 0 and 0,77, and, in addition, it vanishes for 
0- 7 t/ 2, Thus, equation (53) represents a solution for the spherical triangle 
defined by the domain 0 < i9 ^ 0 <2tt; 0 < 0 ^ 77 / 2 . Equivalently stated, this 
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triangle Is bounded by two meridional arcs and the equator, specified by setting 
a ■= b = 0 in inequality (9). In this special case also, the eigenfunction (53) may 
be shown to satisfy the elliptic partial differential equation (12) by direct sub- 
stitution. The boundary conditions (11b) are fulfilled for the particular spherical 
triangle described above. Hence, the eigenfunction (53) is an exact analytic solu- 
tion to problem (11), where 



for any 0 in the interval 0 < 0 < 27t. The original restriction of the solution (51) 
that k = n - 1 = vr/0 be a positive integer may be removed in this special case 
as well. 

Consideration of further special cases obtained when k=n-2,n-3,... 
does not lead to spherical triangles which can be readily identified, so this course 
will not be pursued. However, it is of interest to consider the two special cases 
described above in the event that 0 = 2 t 7 ^, a possibility that was specifically ex- 
cluded in previous discussion. With 0 = 2 t?-, the solution (52) becomes 

u((9, 0) = A s in ^ (sin 

and the corresponding eigenvalue is A. = 3/4. Also, the two-sided spherical 
wedge^'Shaped domain now becomes the entire surface of the sphere except for 
a boundary great circle arc extending between the north and south poles. This 
particular domain may be viewed equivalently as the exterior of a spherical 
triangle in the limit when 0=0. Such an equivalent view will be of significance 
in the discussion of Section VI. In the second special case, the solution (53) 
becomes u(5 ,0) = A sin(0/2)(sin<j5))^^^ cos and the corresponding eigenvalue 
is A. = 15/4. In this instance, the spherical triangle domain degenerates to the 
entire upper hemisphere except for a boundary great circle arc extending from 
the north pole to the equator, and the domain is defined by the values a = b = 0, 

0 = 27T. 



SECTION V 


DISCUSSION OF NUMERICAL RESULTS 

The iterative solution of the finite difference equations by successive over- 
relaxation has been adapted to an electronic digital computer program for cal- 
culating the fundamental eigenvalue of a spherical triangle. The mathematical 
and theoretical basis for the solution has been discussed previously in Sections 
n and in, and the algorithm for the calculation is presented in Appendix A. In 
the present section, numerical results for a number of cases of Interest are 
presented and discussed. 

The first series of spherical triangles to be considered consists of triangles 
containing two right angles where the boundary great circle arcs are formed by 
the equator z = 0 and two meridional arcs. In terms of the spherical triangle 
parameters, such a triangle is specified by a = 0, b = 0 and 0 , where 0 < 0 < 2n. 
It will be recognized that this particular case was solved analytically in Section IV 
for the fundamental eigenvalue 





Table 1 displays calculated values for the fundamental eigenvalue \ for each of 
three values for the spherical triangle parameter © and for various rectangular 
grid parameters, N^j and N^. All the grids utilized a grid spacing parameter 
7=1, i.e., the grid spacing in both the 9 - and ^-directions was constant, as 
specified in equations (13) and (14). The criteria adopted for iterative conver- 
gence for use in this table and in all other numerical results discussed herein 
are ® 10 * • For comparative purposes, the exact value for the 

fundamental eigenvalue, as determined from the analytic relationship \ =\(©) is 
also shown in Table 1. Other columns in Table 1 include the calculated value for 
the optimal scalar relaxation factor cj, the required number m of outer iterations 
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Table 1 

Fundamental Eigenvalues of Double Right Spherical Triangles (a=b=0) 
Using Constant Grid Spacing in the c?5>- Direction (y=l) 


Spherical 

Triangle 

Parameter 

0 

Rectangular 

Grid 

Parameters 

Fundaniental 

Eigenvalue 

Optimal 

Scalar 

Relaxation 

Factor 

O) 

Required 
Number of 
Outer 
Iterations 
m 

Required 
Number of 
Gauss-Seidel 
Inner 
Iterations 

\ 

Required 
Number of 
S.O.R. Inner 
Iterations at 
Convergence 
k 

m 

N, 


Calculated 

Value 

Exact 

Value 

IT 

10 

10 

5.831 

6 

1.460 

8 

85 

24 

77 

20 

20 

5.958 

6 

1.695 

8 

100 

48 

77 

30 

30 

5.981 

6 

1.796 

8 

100 

71 

77 

40 

40 

5.990 

6 

1.860 

8 

100 

99 

77/2 

10 

10 

11.628 

12 

1.506 

8 

100 

27 

77/2 

20 

20 

11.907 

12 

1.723 

10 

100 

49 

77/2 

30 

30 

11.959 

12 

1.807 

11 

100 

73 

77/2 

40 

40 

11.977 

12 

1.850 

11 

150 

97 

77/3 

10 

10 

19.338 

20 

1.518 

12 

100 

27 

77/3 

20 

20 

19.836 

20 

1.729 

13 

100 

50 

77/3 

30 

30 

19.927 

20 

1.806 

13 

100 

78 

77/3 

40 

40 

19.959 

20 

1.846 

13 

150 

110 
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to achieve convergence in the required number kj = K of inner iterations to 
achieve convergence of U by the Gauss-Seidel method for the first outer itera- 
tion, and the required number k^ of inner iterations to achieve convergence of 
U by successive overrelaxation (S.O.R.) for the final outer iteration. The upper 
limits on the number of iterations permitted for inner and outer convergence 
were initially selected to be k^^^^ = 100 and m^^^ = 50, respectively. This value 
for m^a^ was always found to be more than sufficient to promote outer conver- 
gence. However, the value of k^^^^ = 100 was not always sufficient to permit 
inner convergence by successive overrelaxation. (It will be recalled from the 
discussion in Section DI that non-convergence of the Gauss-Seidel iterations 
during the first outer iteration is disregarded.) In such cases of insufficiency, 
the value of k„,a^ was Incremented by 50 successively until inner convergence 
resulted. The number of increments A k^^^^ = 50 necessary to promote inner 
convergence can be ascertained by examination of the column of values for kj. 

In Table 1, it is seen that k^^, = 100 was sufficient for inner convergence except 
for the two = 40 grids when © = 77/2 and 0 = 77 / 3 , in which cases the in- 

crease to kmax ® 150 was then sufficient to promote inner convergence. Further 
examination of this same column of values reveals that only in one instance (for 
the = 10 grid when 0 = 77) did the Gauss-Seidel iterations converge for 

ki < k^^^. In all other cases, k^ = k^^^ , indicating non-convergence of the 
Gauss-Seidel iterations during the first outer iteration. (An alternate, but less 
likely, possibility is that the Gauss-Seidel iterations converged on precisely the 
final permitted iteration, number 

The significant results provided by Table 1 include the qualitative fact that 
as the rectangular grid becomes finer (as measured by an increase in the param- 
eters N |9 and N^), the calculated value for the fundamental eigenvalue \ becomes 
more accurate. Quantitatively speaking, the relative error in \ , given by 



31 


^^xact - Calculated )/Cxact ‘ approximately 0.03 for the = 10 grids, 

0.008 for the = 20 grids, 0.003 for the = 30 grids, and 0.002 for 

the N|j = = 40 grids. These relative errors apply approximately for the three 

values of the spherical triangle parameter 0 , but the relative errors, for a given 
set of grid parameters, are smallest for the case 0 =77 and largest for the case 
0 = 7t/ 3, Of course, use of a finer grid results in a more accurate eigenvalue 
but also entails a substantially greater number of calculations to complete the 
successive overrelaxation method. Note also that the convergence to the exact 
value of X with increasing grid parameters occurs monotonically for each of the 
three 0 values. A similar phenomenon occurs in the calculated values for the 
optimal scalar relaxation factor viz., the values of oj increase with the in- 
creasing grid parameters. Finally, this same phenomenon is repeated in the 
required number of inner and outer iterations. However, the increase in re- 
quired outer iterations m with increasing grid parameters is negligible, while 
the increase in required inner iterations just prior to convergence is con- 
siderable. Note also that the required number of inner iterations for the 
final outer iteration is generally substantially fewer than the required munber 
k j of inner iterations for the first outer iteration. This last characteristic may 
be ascribed to the fact that the convergence criteria are invariant (equal to the 
values eiNNER “ Suter * 10"®) throughout the iterative process. Most of the 
other phenomena delineated above are due to the fact that convergence occurs 
more slowly for a finer grid than for a coarse grid. 

The next series of results, as displayed in Table 2, is intended to investigate 
the effect of adopting a variable grid spacing in the 0 -direction, as discussed 
previously in Section n. For this purpose, a value for the grid spacing parameter 
of 7 - 2 was chosen. This produces a grid in which the density of grid points in- 
creases toward the poles (recall that a singularity occurs at the north pole, 0=0) 



Table 2 

Ftindamental Eigenvalues of Double Right Spherical Triangles (a=b=0) 
Using Variable Grid facing in the 0-Direction {y = 2) 


Spherical 

Triangle 

Parameter 

0 

Rectangular 

Grid 

Parameters 

Fundamental 

Eigenvalue 

Optimal 

Scalar 

Relaxation 

Factor 

Required 
Number of 
Outer 
Iterations 
m 

Required 
Number of 
Gauss -Seidel 
Inner 
Iterations 

Required 
Number of 
S.O.R. Inner 
Iterations at 
Convergence 
k 

m 



Calculated 

Value 

Exact 

Value 

77 

10 

10 

5.750 

6 

1.494 

8 

100 

25 

77 

20 

20 

5.931 

6 

1.716 

8 

100 

49 

77 

30 

30 

5.969 

6 

1.807 

9 

100 

73 

TT 

40 

40 

5.983 

6 

1.852 

9 

150 

97 

n/2 

10 

10 

11.368 

12 

1.519 

10 

100 

26 

7T/2 

20 

20 

11.798 

12 

1.730 

11 

100 

49 

7T/2 

30 

30 

11.910 

12 

1.811 

11 

100 

73 

n/2 

40 

40 

11.949 

12 

1.847 

11 

150 

104 

77/3 

10 

10 

19.123 

20 

1.523 

11 

100 

26 

77/3 

20 

20 

19.609 

20 

1.731 

14 

100 

49 

77/3 

30 

30 

19.825 

20 

1.801 

13 

100 

81 

77/3 

40 

40 

1 

19.902 

20 

1.847 

13 

; 150 

105 
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and decreases toward the equator, in a fashion symmetrically arranged about the 
equator. The same spherical triangles and the same values for the rectangular 
grid parameters are considered as in Table 1. The outcome of using this par- 
ticular variable grid spacing on the calculated values of X. is seen to be larger 
relative errors in each of the 12 cases considered, as compared to the results 
of Table 1 for a constant grid spacing. Hence, this method of varjdng the grid 
spacing in the 0 -direction in order to account for the singularity in 4> does not 
achieve the desired effect. Virtually all other qualitative characteristics noted 
in the Table 1 results also apply to Table 2, however. In particular, more re- 
fined grids produce more accurate calculated values for the fundamental eigen- 
value, and the convergence to the known exact value for ^ occurs monotonically 
in each of the three 0 cases. 

A variable grid spacing in the <?b -direction based upon the parametric value 
7 = 1/2 was next investigated, and the results are displayed in Table 3. This 
value of y produces a grid in which the density of grid points increases toward 
the equator and decreases toward the poles, again in a manner symmetrical with 
respect to the equator. The same spherical triangles with two right angles and 
the same values for the rectangular grid parameters are considered as in Tables 
1 and 2. This particular variable grid spacing produces calculated values of ^ 
with larger relative errors than both Tables 1 and 2 for the 0 = 77 cases and for 
the 0 = 7t/2 cases for the coarser = 10 and = 20 grids. How- 
ever, for the finer = 30 and = 40 grids when 0 = 77 / 2 , the cal- 

culated X. are more accurate in the y = 1/2 grid spacing (Table 3) than in the 
7 = 2 grid spacing (Table 2), though not so accurate as in the constant grid spac- 
ing (Table 1). For the case 0 = tt/3, the result for k in Table 3 is least accurate 
for the coarsest = 10 grid, but most accurate for the finer three grids, 

even when compared against the constant grid spacing results. For 0 = tt/3, 



Table 3 

Fundamental Eigenvalues of Double Right Spherical Triangles (a=b=0) 
Using Variable Grid Spacing in the ^-Direction (7 = 1/2) 


Spherical 

Triangle 

Parameter 

0 

Rectangular 

Grid 

Parameters 

Fundamental 

Eigenvalue 

K 

C^timal 

Scalar 

Relaxation 

Factor 

CO 

Required 
Number of 
Outer 
Iterations 
m 

Required 
Number of 
Gauss-Seidel 
Inner 
Iterations 

kj 

Required 
Number of 
s.o.B. Iimer 
Iterations at 
Convergence 
k 

tn 

N, 


Calculated 

Value 

Exact 

Value 

7 T 

10 

10 

5.118 

6 

1.427 

8 

73 

22 

77 

20 

20 

5.592 

6 

1.673 

8 

100 

45 

V 

30 

30 

5.773 

6 

1.780 

8 

100 

62 

7 T 

40 

40 

5.856 

6 

1.844 

8 

100 

85 

77/2 

10 

10 

10.692 

12 

1.488 

10 

98 

26 

77/2 

20 

20 

11.739 

12 

1.713 

11 

100 

46 

77/2 

30 

30 

11.919 

12 

1.797 

11 

100 

75 

77/2 

40 

40 

11.969 

12 

1.850 

11 

100 

94 

77/3 

10 

10 

19.061 

20 

1.508 

11 

100 

27 

77/3 

20 

20 

19.888 

20 

1.720 

13 

100 

52 

77/3 

30 

30 

19.961 

20 

1.806 

13 

100 

74 

77/3 

40 

40 

19.979 

20 

1.848 

13 

150 

104 
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the relative error in \ is approximately 0,047 for the * 10 grid, 0.006 

for the Ng = = 20 grid, 0.002 for the = 30 grid, and 0.001 for the 

= 40 grid. These relative errors in k are considerably smaller than 
those for the corresponding grids in the cases © = tt/2 and 0 = 77 , The conclu- 
sion appears to be that variable grid spacing of the type used in the Table 3 
results can produce a more accurate calculated value for the fundamental eigen- 
value than that produced by use of constant grid spacing, but the relative accu- 
racies depend upon the characteristics of the spherical triangle under consider- 
ation. Other convergence qualities noted in the previous Tables 1 and 2 results 
for the remaining parameters displayed also hold for the Table 3 results. 

The abbreviated display in Table 4 considers the special case 0 = 2 t 7 - only, 
in which the spherical triangle degenerates to the upper hemisphere with bound- 
aries consisting of the equator and one meridional great circle arc (a ’’slit") 
extending from the north pole to the equator. This degeneracy was specifically 
considered at the end of Section IV, and the analytic relationship 



yields an exact value for the fundamental eigenvalue of 15/4. The same sequence 
of values for the rectangular grid parameters is adopted in Table 4 as in the 
earlier tables, and constant grid spacing in the 0 -direction with 7 = 1 is utilized, 
as was the case for Table 1 . Once again, the finer grids are seen to produce 
more accurate calculated values of k with convergence to the known exact value 
occurring monotonically, but in this instance the convergence is monotone de- 
creasing. In the previous results of Tables 1 , 2, and 3 , the convergence was 
always monotonic and increasit^. Otherwise, the results of Table 4 show the 
usual convergence properties noted earlier and serve pr ima rily to verify the 
analytic solution for a particular special case. 



Table 4 

Fundamental Eigenvalue of Slit Upper Hemisphere (a=b=0; 0 = 2tt) 
Using Constant Grid facing in the qfc -Direction (7 = 1) 


Spherical 

Triangle 

Parameter 

0 

Rectai^lar 

Grid 

Parameters 

Fundamental 

Eigenvalue 

Optimal 

Scalar 

Relaxation 

Factor 

CO 

Required 
Number of 
Outer 
Iterations 
m 

Required 
Number of 
Gauss-Seidel 
Inner 
Iterations 
kx 

Required 
Number of 
S.O.R. Inner 
Iterations at 
Convergence 
k 

m 

N, 


Calculated 

Value 

Exact 

Value 

2t7 

10 

10 

3.845 

3.75 

1.385 

10 

59 

21 

2t7 

20 

20 

3.826 

3.75 

1.647 

10 

100 

42 

2tt 

30 

30 

3.806 

3.75 

1.760 

10 

100 

59 

2tt 

40 

40 

3.794 

3.75 

1.837 

9 

100 

83 
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The next series of results in Table 5 investigates another of the special cases 
discussed previously in Section IV, that of the spherical wedge. It will be recalled 
that the spherical wedge is a degenerate case of the spherical triangle which is 
bounded by two meridional great circle arcs extending from pole to pole, and it is 
specified by the limits a = 0 , b ■* - when 0 < tt . The exact value for the funda- 
mental eigenvalue is given by 



When 9. = 77 with a= 0,b ^ the spherical wedge becomes a hemisphere, with 
fundamental eigenvalue A. = 2. However, when 0 assumes a value in the range 
TT < 0 < 277 with a = 0, b - -00, the domain specified is still the hemispherical 
wedge, with the corresponding eigenvalue A. = 2. As can be seen from inequality 
(9), the only value of 4> in the domain for a = 0, b ^ - oo^ and d in the range 
7T< 6 < 2tt is 4> => 0, the north pole. For numerical purposes, the value adopted 
for b in the Table 5 results was -10® . The standard sequence of values for the 
rectangular grid parameters is adopted as in the earlier tables, and constant 
grid spacii^ is utilized. For both the quarter-spherical wedge {0 » 77 / 2 ) and 
the hemispherical wedge (0 = 77), the convergence to the known exact value of \ 
is monotonically increasing, as in the first three tables, with the more refined 
grids. Also, the calculated values for the fundamental eigenvalue are seen to 
possess considerable accuracy for all grids. Now in the cases 0=3 n/2 and 
0 = 2 t7 , the convergence in A. is also monotonically increasing with the refine- 
ment in the grids, but the corresponding relative errors are larger than for the 
smaller values of 0 . This is undoubtedly due to the fact that a larger propor- 
tion of the grid points in the 0 = 877 / 2 , 2 77 cases are "superfluoi^," i.e., propor- 
tionately fewer grid points occur interior to the spherical wedge domain than in 
the corresponding cases for 0 = 77/2 and 0 = 77 , for a given set of rectangular 



Table 5 

Fundamental Eigenvalues of Spherical Wedges (a=0, b- - <») 
Using Constant Grid Spacing in the ^ -Direction (y = 1) 


Spherical 

Triangle 

Parameter 

0 

Rectangular 

Grid 

Parameters 

Fundamental 

Eigenvalue 

Optimal 

Scalar 

Relaxation 

Factor 

CtJ 

Required 
Number of 
Outer 
Iterations 
m 

Required 
Number of 
Gauss -Seidel 
Inner 
Iterations 

Required 
Number of 
S.O.R. Inner 
Iterations at 
Convergence 
k 

m 



Calculated 

Value 

Exact 

Value 

■n/2 

10 

10 

5.926 

6 

1.543 

8 

100 

29 

77/2 

20 

20 

5.982 

6 

1.750 

8 

100 

54 

77/2 

30 

30 

5.992 

6 

1.808 

9 

100 

95 

77/2 

40 

40 

5.995 

6 

1.854 

9 

150 

122 

7J 

10 

10 

1.977 

2 

1.572 

6 

100 

31 

U 

20 

20 

1.994 

2 

1.774 

6 

100 

59 

77 

30 

30 

1.997 

2 

1.824 

6 

150 

101 

77 

40 

40 

1.999 

2 

1.838 

6 

200 

177 

377/2 

10 

10 

1.828 

2 

1.512 

6 

100 

26 

3t7/2 

20 

20 

1.852 

2 

1.722 

6 

100 

48 

3t7/2 

30 

30 

1.917 

2 

1.830 

6 

100 

73 

377/ 2 

40 

40 

1.961 

2 

1.848 

6 

150 

90 

2t7 

10 

10 

1.929 

2 

1 

1.465 

6 

83 

23 

2tt 

20 

20 

1.972 

2 

1.686 

6 

100 

45 

2tt 

30 

30 

1.982 

2 

1.788 

6 

100 

62 

2jt 

40 

40 

1.987 

2 

1.845 

6 

150 

83 


CO 

00 
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grid parameters. Nonetheless, convergence to the proper exact value of K does 
occur monotonically for all 0 values in Table 5. Convergence properties of 
other parameters displayed in the table are similar to those observed in pre- 
vious results. 

This concludes the discussion of numerical results obtained for all the 
special cases considered in Section IV for which analytical solutions were ob- 
tained, thereby providing known exact values for the fundamental eigenvalues. 

A final series of results is presented in Table 6 for an archetypical case ( a so- 
called oblique spherical triangle) representing a numerical solution to the gen- 
eral problem for which the exact fundamental eigenvalue is not known. A family 
of spherical triangles defined by the parameters a = 0, b = -1 is considered, 
where 0 assumes each of the four values 7 t/2, n, 377 / 2 , and 2 t 7. The oblique 
plane defining the non-meridional side of the spherical triangles, as given by 
equation {7c), is z = -y, and this plane forms an angle of 77/4 with the equator. 

The great circle determined by this oblique plane intersects the equator at 
0, 77, and 2 77, and the interior domains of the spherical triangles appear as 
shown in Figure 5. As in the results discussed previously, a standard set of 
values for the rectangular grid parameters and constant grid spacing defined 
7=1 are utilized. In three of the four © cases displayed in Table 6, the con- 
vergence of the calculated values of X to an unknown exact fundamental eigen- 
value is monotonically increasing with the refinement in grids, but in the anom- 
alous 0 = 77 case, oscillatory behavior in the calculated values of \ is observed. 
No explanation for this lack of monotonicity is offered here, but this phenomenon 
has been observed elsewhere (Reference 6, pp. 351-352) in the numerical calcu- 
lation of fundamental eigenvalues. Furthermore, in the cases 0 =877/2 and 
0 = 277 , although the convergence shown is monotonic, it certainly is not uniform 
to the extent generally seen in earlier results for the known solution special cases. 



Table 6 

Fundamental Eigenvalues of Oblique Spherical Triangles (a=0,b=-l) 
Using Constant Grid Spacing in the <?b-Direction (y = 1) 


Spherical 

Triangle 

Parameter 

Rectangular 

Grid 

Parameters 

Calculated 

Fundamental 

Optimal 

Scalar 

Relaxation 

Required 
Number of 
Outer 

Required 
Number of 
Gauss -Seidel 
Inner 

Required 
Number of 
S.O.R. Inner 
Iterations at 

0 

N, 


Eigenvalue 

X 

X* actor 

Ct) 

It 0 r ati OHS 
m 

Iterations 

Convergence 

k 

m 

'/t/2 

10 

10 

7.392 

1.520 

8 

100 

28 

7T/2 

20 

20 

7.569 

1.731 

10 

100 

49 

n/2 

30 

30 

7.633 

1.809 

11 

100 

77 

n/2 

40 

40 

7.668 

1.854 

11 

150 

100 

77 

10 

10 

2.795 

1.519 

6 

100 

28 

77 

20 

20 

3.027 

1.724 

7 

100 

52 

77 

30 

30 

3.000 

1.823 

8 

100 

76 

77 

40 

40 

3.061 

1.843 

8 

150 

118 

377/2 

10 

10 

2.435 

1.476 ; 

8 

90 

24 

3tt/2 

20 

20 

2.615 

1.697 

7 

100 

47 

377/2 

30 

30 

2.631 

1.796 

8 

100 

66 

377/2 

40 

40 

2.664 

1.844 

8 

150 

88 

2t7 

10 

10 

2.396 

1.435 

8 

73 

22 

277 

20 

20 

2.581 

1.665 

8 

100 

42 

2t7 

30 

30 

2.597 

1.771 

8 

100 

57 

2t7 

40 

40 

2.652 

1.811 

9 

100 

92 







Figure 5. Interior domains of oblique spherical triar^les defined by the 
parametric values a = 0, b= -1; 0 = tt/2, tt , 37t/2, and 2t7 (indicated by 
the shaded regions to the left of each of the four respective ^ = 0 values) 


>(>• 

|m^ 
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The non-uniform convergence of Table 6 may well be due to the fact that the 
oblique spherical triangles contain a boundary great circle arc which appears 
as a transcendental curve in the two-dimensional 6 ,0-plane (see Figure 5). In 
the general problem, this transcendental curve will be present, although it did not 
appear in the earlier results of this section for the special cases when all spher- 
ical surface boundaries coincided with constant- value 6 and 0 great circle arcs 
(i.e., lines of spherical longitude and latitude). The presence of a curved bound- 
ary in the 8 ,0 -plane requires that a series of straight line segments connecting 
appropriate grid points <9. , 0^ be utilized to approximate the southern boundary 
of the spherical triangle. As the grid parameters are changed, the pre- 

cise location of the approximate triangle boundary varies. This could perhaps 
explain the non-uniform convergence in the calculated values of \ . The conver- 
gence characteristics of the relaxation factor and the required numbers of 
iterations as shown in Table 6 are similar to those noted in the earlier results 
for the special cases. 



SECTION VI 


PROBLEM OF A SPHERICAL SURFACE WITH A SLIT 
A« The Slit Domain 

In this section, the problem of determining the fundamental eigenvalue for a 
domain which is the exterior of a spherical "triangle" in the degenerate case of 
® = 0 will be considered. In such a limiting case, the triangle degenerates to a 
meridional great circle arc or "slit." This problem is the second of the two 
cases mentioned in the introductory section for which the second-order elliptic 
partial differential equation (12) will be solved by approximation methods, with 
the boundary conditions supplied by equation (lib). 

The slit domain D may be specified by a single slit extension parameter 
which shall be denoted $ , where 0 < # < tt. The slit may be readily placed in a 
canonical position so that it extends from the spherical north pole (intersection of 
the z-axis with the unit sphere, as previously) southward along the arc of a great 
circle. In terms of the spherical triangle parameters utilized previously, the 
value for $ may be determined from equation (8) by setting 0 = 0 (hence 6 - 0 ) 
so that the resulting value is equal to It is then seen that $ = arccot a = 
arcsin (a^ + 1)'^''^ , with b arbitrary. 

The spherical surface may be mapped directly onto the 0,0 -plane in the 
usual manner, as shown in Figure 6. Note that the domain of interest now in- 
cludes the entire spherical surface, 0 ^ 4 >^ tt \ 0 ^ 9 ^ 2t 7. A network of grid 
points <9. , , as defined by equation (13), is superposed over the spherical sur- 

face in the manner used previously, but now the relations between the number of 
grid intervals and the constant grid spacings are 



SPHERICAL 

SLIT 



Figure 6. Network of grid points , <p. superposed over entire spherical 
surface with a slit mapped onto a plane rectangular region 
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In the 0 , c^-plane, the single slit appears along some portion of three of the rec- 
tangular boundaries, viz., the entire line representing the north pole = 0, a 
portion of the line representing the slit meridian (9 = 0, and a like portion of the 
line i9 s! 2 t 7 also representing the same slit meridian. Note that, in general, the 
slit boundary will not terminate on a grid point. 

B. Boundary Conditions 

The boundary conditions (11b) to be imposed on the discretized version of 
the eigenvalue equation for the problem of the spherical surface with a slit may 
be deduced from Figure 6. They are 

^ 2, N^, 

and 

^o.i = \,j = 0 for j = 0. 1, 2,.... N^, if(?^. <$. (56) 

The boundary conditions (55) refer to the north pole of the sphere, and the condi- 
tions (56) refer to grid points coinciding with the slit meridian boundary. Note 
that, in general, there are non-zero values of the grid function Uj . occurring 
along the rectangular grid boundary in Figure 6; in particular, such values may 
appear along the slit meridians beyond the extension of the slit (i.e,, U and 

0, i 

Uj^ , for which 0. > $) and at the south pole (i.e., U. for i * 0, 1, 2 N ). 

In contrast, all grid function values occurring along the rectangular grid boundary 
in the problem of the spherical triangle vanished. The finite difference equation 
(22) , derived as an approximation to the continuous partial differential equation 
(12) for the problem of the spherical triangle, remains valid for all rectangular 

interior grid points ( 9 ., 0. ), i = 1, 2 N^ - 1; j = 1, 2 N^ - 1 in the 

problem of the spherical surface with a slit. However, it must be modified in 
order to apply- to the rectangular boundary grid points represented by the sub- 
script values i = 0, N^ ; j = 0, 1, 2 N^ if 0. > # and i = 0, 1, 2, . . . ,N^; 

j = N^. 
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Since the rectangular boundary grid points along the line i = 0 coincide with 
the respective rectangular boundary grid points aloi^ the line i = insofar as 
the spherical surface geometry is concerned, it is evident that the following 
periodic boundary conditions must be imposed in addition to conditions (55) and 
(56): 

j =0- 2. --. N^. (57) 

In view of this periodicity, the finite difference equation (22) may be extended to 
the line i = by the simple device of replacing the undefined grid components 
Un^+ 1, j I j = 0, 1, 2, ...» by the respective defined grid components j , 

j = 0, 1, 2 N^, whenever the former occur. That is, for i = N^, equation 

(22) becomes in its modified form: 




\e.V^ for j = 1, 2, - 1. if 0. > 


(58) 


where a . , b , , c . , and e . are defined as previously. Now the finite difference 
equation (22) may be further extended to the line i = 0 merely by invoking equa- 
tion (57), In order to complete the extension of the finite difference equation to 
the rectangular boundary grid points, it is necessary to derive relations for the 

grid points, U. , i = 0, 1, 2, . . . , N5, at the south pole. 

* ^ ^ 

C. Difference Equations at the South Pole 
Difference equations pertaining to the grid points at the south pole may be 
derived by integrating the partial differential equation (12) over the disc O<0 <2tt; 

^ ‘Pm - 7T(see Figure 7). Note that the first term on the left in equa- 
tion (12) vanishes when integrated over the disc because of the 2 n periodicity in 


the & variable. After integrating the second term on the left in equation (12) 
explicitly with respect to the variable, it is seen that 
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^3N^/4 =37T/2 

Figure 7. Disc (indicated by shading) utilized to derive difference 
equations valid at the spherical south pole 
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u s in (pdcpdd ~ 0- 


^N^-l/2 


(59) 


The first term on the left in equation (59) may be approximated, by use of the cen- 
tered first difference quotient (19), as a finite sum of the form, 


N, 


-z: 




h sin ch 


i = 1 


^jb.N- 


Na- 1/2 ' 


(60) 


The second term on the left in equation (59) may be approximated by 


In equation (61) and equations to follow, it is implicitly understood that U. 


(61) 


i.N. 


represents the unique approximate value of the solution u at the south pole for 
all values of the subscript i. That is, 


= i = 0- 1- 2,..., N^. 


(62) 


Now if the approximations (60) and (61) are utilized in equation (59) and the re- 
sulting equation multiplied by the factor (2/hg ), it is seen that 


2 s i n 0^, -1 /„ 


%-l/2 




i = l 


Define the hitherto undefined quantity e by 


4t7 

- (1 + COS d> V 

hg ' 1/2 


( 64 ) 


which, by use of equations (54a) and (17), may be rewritten as 
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Eciuation (63) may now be written, with the aid of e^atlons (33b) and (64) , as a 
finite difference equation for the grid points U. at the south pole in the form 




( 66 ) 


\ i = i / 

Note that the difference equation (66) is not a five-point equation, as are the dif- 
ference equations (22) and (58), but is instead a (N^ + l)-polnt equation. This 
completes the derivation of the additional finite difference equations and bound- 
ary conditions necessary to extend the solution of the problem of a spherical 
surface with a slit to the complete set of rectangular Interior and boundary grid 
points. 


D, Application of Successive Overrelaxation 
In the matrix formulation of this problem, U represents the ordered set of 
unknown values Uj . written as a vector defined on the network of rectangular 
interior and boundary grid points , 4>-^. The number of components in the vec- 
tor U is equal to the number of rectangular interior grid points, (N. - 1)(N^ - 1), 
plus a number of rectangular boundary grid points equal to the number of grid 
points along either of the lines i » 0 or i = Ng which do not coincide with the 
slit, including the single point at the south pole. This number of components of 
U is equal to the number of independent unknowns U. . to be determined from 
the three finite difference equations (22), (58), and (66). Equivalently, the num- 
ber of components of U is equal to the total number (N^ + 1)(N^ + 1) of rectangu- 
lar interior and boundary grid points reduced by the constraints of the boundary 
conditions at the north pole (55), along the slit (56) and its meridional extension 
(57), and at the south pole (62). Then the matrix form (29) of the finite difference 
equations, with the boundary conditions included, still applies, where A and E 
are now square matrices of an order equal to the dimensionality of U and with 
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somewhat different properties than those enumerated in the problem of the spher- 
ical triangle. 

An example of the matrix formulation for the special case of = 4, = 6 

will now be considered, where $ = tt/ 4 and 7=1 are assumed, so that constant 
grid spacing in the 4 > -direction results. Note that this value of $ results in a 
slit boundary whose termination does not coincide with a grid point. Figure 8 
displays the network of grid points, the location of the slit boundary, and a natural 
ordering of the grid points in which successive horizontal grid lines are scaimed 
in order of increasing i index with the j index increasing on each new line. The 
rectangular boundary grid points are chosen along the line i = and the single 
point at the south pole is arbitrarily placed at the position (1, N^). The vector 
U in this case has dimensionality 20, of which 15 components are associated 
with rectangular interior grid points and the 5 remaining components with rec- 
tangular boundary grid points. Figure 9 provides the square matrix A of order 
20, the square diagonal matrix E of order 20, and the vector U for this special 
case. 

In general, the matrix A is sparse, with positive diagonal entries and non- 
positive off-diagonal entries. Because of the fact that the three finite difference 
equations (22), (58), and (66) are symmetric, the matrix A is also symmetric. 
However, A no longer retains the block tridiagonal property that it previously 
exhibited in the problem of the spherical triangle. The matrix E is diagonal 
with positive entries, by virtue of equations (24c) and (65). The elements of both 
matrices A and E depend only on the j index and are independent of i. 

As an Initial estimate for the power method iterations (34), is defined 

so that all components are assigned the value unity, except for the zero compo- 
nents associated with grid points at the north pole and along the slit boundary, 
in accordance with the boundary conditions (55) and (56). In particular, components 
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SOUTH POLE 

(4,6) 


(4,5) 


(4,4) 


(4,3) 


(4,2) 


(4,1) 


(0,0) (1,0) (2,0) (3,0) (4,0) 

0 = 0 0=iL d^ir 0 = 0 = 2ir 

2 2 

Open circles represent grid points along slit boundary (heavy line) corresponding to zero values of the grid function U. 

Filled circles represent grid points corresponding to non-zero values of U. 

Numbers in parentheses indicate grid point (i, j) co-ordinates. 

Numbers In circles indicate natural ordering of the grid points. 
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Figure 8. Network of grid points for the special case 
= 4, N^ = 6, $ = 7 t/ 4, and7= 1 
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Figure 9. The square matrix A of order 20, the vector U of 20 components 
(the superscript indicates transpose), and the square diagonal matrix 
E of order 20 for the special case = 4, 6, $= tt/ 4, and 7 » 1 
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associated with the grid points at the south pole and along portions of the slit 
meridian beyond the extension of the slit are initialized to the value unity. 

The iterative method of point successive overrelaxation as applied to the 
finite difference equation (22) produces iterates whose components asso- 

ciated with rectangular interior grid points are given by equation (36). For values 
of the j index such that 0. > f , the iterates U whose components are associ- 
ated with rectangular boundary grid points along the line i = are given instead 
by 


= - [b. + c. + b. . 

N^ij g jlij j 

+ c : .].+ (1 

J “ 1 . j - 1 Na , J ' Nj9 , J 


for j ^ 1, 2, - 1, if 0. > $; k = 1, 2, 


(67) 


The iterate components in equation (67) are, of course, based upon the finite dif- 
ference equation (58). For j = , the iterates whose component is asso- 

ciated with the grid points at the south pole are given by 


U(k+i> _ ^ 


i ,N 




2L, 

i= 1 


i.Nu 




+ (1 1 , 2 ,... 


These iterates are based upon the finite difference equation (66), and because of 
the uniqueness specified in equation (62), it may be seen that equation (68) is 
actually independent of the i index. Singularity problems do not arise, since 
Ng, and c never vanish, by definition and equation (24b), respectively. 
Furthermore, it is important that the iterate components given by equations (36), 
(67), and (68) be computed according to the natural ordering of the grid points 
associated with the components of the vector U- j . Use of the natural ordering 
will insure that the proper values of the superscripts k and (k+ 1), specifying 
the inner iteration number, occur automatically in the calculation process. 
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The vectors F that appear as k-independent constants in the inner iteration 
equations (67) and (68) are defined, in accordance with the finite difference equa- 
tions (58) and (66), merely by extending the definition (37) in the natural way to 
include components associated with rectangular boundary grid points. That is, 
with use of definition (65) , 


F. 


i.j " 2 N^; j = 1. 2,.... N,. 


(69) 


Similarly, the definition (38a) of the convergence parameters must be 

extended to the full range of grid point indexing as follows: 


e(‘‘) = 


max 


uCk+i) k= 1, 2... 

1 , j 1 * J I * ‘*■1 


r(k) 


(70) 


The range of the summations involved in the Rayleigh quotient (35) must likewise 
be adjusted to account for the proper range of grid point indexing. The vector 
inner products may be written in component form as - P/Q, where, in con- 
trast to equations (42) , 


and 


P - (U<"’>, EUC 




i=l \ j= 1 



(71a) 


Q = (A'iEU("'). EU("’>) 


z: fc 


F: 






0 


(71b) 


i=i \i = i / 

In the summations appearing in equations (71), components associated with rec- 
tangular boundary grid points along the line i * 0 are excluded, since such com- 
ponents are redundant to those along the boundary line i » , the latter of which 

are included in the summations. Also, the single components associated with the 
grid points at the south pole are presented separately as the final term on the 
right side of equations (71). Finally, the component form (45) of the power method 
must be extended to the full range of grid point indexing as in the following: 



55 


i. j 


for i = 0, 1, 2,..., N^; j = 1, 2, . . . , 


(72) 


The theoretical demonstration of convergence of the successive overrelaxa- 
tion iterates (36), (67), and (68) for any initial vector selected follows from 

the Ostrowski-Reich theorem (Reference 7, p, 123) for a relaxation factor w in 
the range 0 < u> <2. As stated previously in Section HI, this theorem can be in- 
voked because A may be shown to be a symmetric irreducibly diagonally domi- 
nant matrix which is, therefore, positive definite. Although the method of suc- 
cessive overrelaxation will converge for all relaxation factors in the above range, 
the most rapid convergence occurs for some optimal value oi of a> in the 

opt 

interval 1 < w < 2. However, since in the problem of a spherical surface with a 
slit, difference equations arise which are not of the five-point type, it is not true 
that the Jacobi matrix associated with A is cyclic of index 2, as was the case 
previously in the problem of the spherical triangle. This means that the theo- 
retical expression (40) for no longer applies. An expression for the actual 
optimal value of w Is not known for the current problem. However, it has 
been suggested (Reference 6, p. 262) that the value of oj produced by the formula 
(40) may be reasonably close to the unknown optimal value , even when the 
theoretical expression (40) does not rigorously apply. For this reason and in 
light of the fact that a better method for approximating the true value of co is 

Op t 

not available, the formula (40), based upon the convergence (39) of parameters 
r defined in equation (38b), is retained for approximating in the nu- 
merical results that follow. 

The implementation of the numerical technique, as previously described in 
Section HI, for the solution of the finite difference equations by point successive 
overrelaxation and the subsequent iterative determination of the fundamental 
eigenvalue proceeds in precisely the same manner for the digital computer 
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solution of the problem of a spherical surface with a slit. Application of these 
techniques is clarified by the algorithmic presentation within Appendix B of the 
methods discussed in the current section. In the algorithm of Appendix B, the 
vector U is considered to be an array with (N^ + 1)(N^ + 1) entries, so that an 
entry of U is associated with each rectangular interior and boundary grid point 
regardless of the value of the slit extension parameter $ . From a computational 
standpoint, this device of enlarging the dimensionality of U, with the constraints 
(55), (56), (57), and (62) imposed upon the entries of the array, is merely a con- 
venience in indexing. 


E. Numerical Results 

In the numerical applications to be described, the slit extension parameter $ 
was assigned a range of values in the interval 0 < $ < tt in a systematic manner, 
and for each such value of $ , several rectangular grid parametric values were 
considered. Table 7 displays calculated values for the fundamental eigenvalue k 
for each of the values $ = n( tt/ 8 ), n * 1, 2 , . . . ,7 and for the familiar gamut of 
rectangular grid parameters, = N^= 10, 20, 30, 40. All the grids utilized a 
grid spacing parameter 7 = 1 , so that uniform 0 ,<p -plane grid spacing resulted. 
The criteria for iterative convergence remain » ^outer ” The 

description provided in Section V of the significance of the various columns within 
the tables applies equally well to Table 7. The value for the increment 
used to promote inner convergence when the current value of k^^^ is not suffi- 
cient remains at 50, as previously. 

The significant results provided by Table 7 include the qualitative fact that 
the calculated value of k increases with the extent of the spherical slit as meas- 
ured by $. However, the convergence in k with increasing rectangular grid 
parameters for a given $ is by no means monotonic, as generally was the case 
In the problem of the spherical triangle. The calculated values for the 



Table 7 

Fundamental Eigenvalues of a Spherical Surface with a Slit 
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Spherical 

Slit 

Extension 

Parameter 

Rectangular 

Grid 

Parameters 

Calculated 

Fundamental 

Eigenvalue 

\ 

Near- 

Optimal 

Scalar 

Relaxation 

Factor 

(a) 

Required 
Number of 
Outer 
Iterations 
m 

Required 
Number of 
Gauss -Seidel 
Inner 
Iterations 

Required 
Number of 
S.O.R, Inner 
Iterations at 
Convergence 
k 

m 

Ns 


77/8 

10 

10 

0.209 

1.794 

4 

100 

90 

77/8 

20 

20 

0.193 

1,900 

4 

200 

195 

77/8 

30 

30 

0.188 

1.937 

5 

300 

283 

77/8 

40 

40 

0.199 

1.960 

5 

400 

317 

77 /4 

10 

10 

0.263 

1.773 

4 

100 

80 

77 /4 

20 

20 

0.273 

1.883 

5 

200 

164 

77/4 

30 

30 

0.259 

1.926 

5 

250 

245 

77 /4 

40 

40 

0.265 

1.953 

5 

350 

271 

877/8 

10 

10 

0.317 

1.755 

5 

100 

73 

377/8 

20 

20 

0.325 

1.874 

5 

200 

150 

377/8 

30 

30 

0.327 

1.917 

5 

250 

223 

377/8 

40 

40 

0.329 

1.949 

5 

300 

242 

7 t /2 

10 

10 

0.427 

1.722 

5 

100 

61 

■n /2 

20 

20 

0.406 

1.860 

5 

150 

133 

77/2 

30 

30 

0.399 

1.910 

5 

200 

198 

7 t /2 

40 

40 

0.396 

1.941 

5 

300 

229 

577/8 

10 

10 

0.486 

1.705 

5 

100 

56 

577/8 

20 

20 

0.465 

1.851 

5 

150 

122 

577/8 

30 

30 

0.457 

1.905 


200 

180 

577/8 

40 

40 

0.469 

1.936 

5 

300 

203 

377/4 

10 

10 

0.550 

1.688 

5 

100 

61 

377/4 

20 

20 

0.563 

1.835 

5 

150 

104 

377/4 

30 

30 

0.544 

1.884 

5 

200 

182 

377/4 

40 

40 

0.551 

1.897 

6 

300 

296 

777/8 

10 

10 

0.619 

1.668 

6 

100 

46 

777/8 

20 

20 

0.637 

1,811 

6 

100 

100 

777/8 

30 

30 

0.644 

1.860 

6 

200 

169 

777/8 

40 

40 

1 0.647 

1.880 

6 

300 

262 
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’’near-optimal** scalar relaxation factor w do increase monotonically with the re- 
finement in the grid, and it is also observed that for a fixed pair of rectangular 
grid parameters , the calculated ct) increases progressively as § decreases. A 
similar monotonic increase in the required number of inner and outer iterations 
occurs with the grid refinement for each value of $. Note also that for a fixed 
pair of Ng, values, the required number of inner iterations at convergence 
increases dramatically as $ decreases* A parallel phenomenon is seen for the 
required number kj of Gauss-Seidel inner iterations for the first outer iterations. 
The change in the required number m of outer iterations, although a much more 
moderate variance, exhibits a reverse trend, viz., a slight decrease with the de- 
crease in $ . In sum, these results indicate that convergence in K occurs more 
slowly (1) for a finer grid than for a coarse grid, and (2) as the slit extension 
decreases toward the north pole. Also, in comparison with the results for the 
problem of the spherical triangle, the Table 7 results indicate that convergence 
in A- to the 10“® criteria requires a substantially greater number of inner 
iterations and a moderately smaller number of outer iterations. The slower 
convergence of the successive overrelaxation method for the problem of a 
spherical surface with a slit is no doubt due to the fact that an expression for 
the optimal value of is not available, as explained previously. 

A further investigation was conducted to determine the effects of utilizing 
values of near the polar limits of 0 and n for the slit boundary. The results 
appear in Table 8. It was desirable to choose two values, # ^ and $ 2 » such that 
0 <• < 7 t/ 40 and tt- {tt/40) < $2 < ^ in order that none of the grid points, 

except those at the north pole, would coincide with the slit boimdary (in the case 
of $ j) and in order that all of the grid points along the slit meridian, except those 
at the south pole, would coincide with the slit boimdary (in the case of § 2 )» for all 
values of the rectangular grid parameters, = N^= 10, 20, 30, 40. The values 
arbitrarily adopted that meet these criteria were $ j = 77/100 and ^*2 ~ 99^/100. 



Table 8 

Fandamental Eigenvalues of a ^herical Surface with a Slit 
Near Slit Boundary Polar LimitiDg Values 


Spherical 

Slit 

Extension 

Parameter 

$ 

Rectangular 

Grid 

Parameters 

Calculated 

Fundamental 

Eigenvalue 

Near- 

Optimal 

Scalar 

Relaxation 

Factor 

CO 

Required 
Number of 
Outer 
Iterations 
m 

Required 
Number of 
Gauss -Seidel 
Inner 
Iterations 

Required 
Number of 
S.O.R, Inner 
Iterations at 
Convergence 
k 

m 

N, 


77/100 

10 

10 

0.147 

1.825 

4 

150 

110 

77/100 

20 

20 

0.123 

1.919 

4 

250 

249 

77/100 

30 

30 

0.112 

1.950 

4 

400 

383 

77/100 

40 

40 

0.105 

1.971 

4 

500 

459 

9977/100 

10 

10 

0.697 

1.646 

6 

100 

40 

9977/100 

20 

20 

0.724 

1.790 

6 

100 

87 

9977/100 

30 

30 

0.733 

1.843 

6 

200 

143 

9977/100 

40 

40 

0.737 

1.863 

6 

300 

222 
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The calculated fundamental eigenvalues display a monotonic decrease for # = 77/IOO 
and a monotonic increase for $ = 99^^/100, with the refinement in grids. More 
significantly, the calculated A. values are in accordance with the trend noted in 
the Table 7 results, i.e., k increases with $. In addition, all the characteristics 
noted in the results of Table 7 for the scalar relaxation factor cd and the required 
numbers of outer and inner iterations m, and are exhibited further by the 
results of Table 8. In particular, the slow convergence of the inner iterations 
for $ * 7t/100 and the correspondingly rapid convergence of the inner iterations 
for $ = 99 77/100 are emphasized in the Table 8 results. 

Finally, the value of the calculated fundamental eigenvalue as a function of 
the spherical slit extension parameter is graphically displayed for three of the 
rectangular grids in Figure 10 {the case = 30 is omitted from the figure 

for purposes of legibility only). It is seen that X versus $ is a very nearly 
linear relationship, with a slope d^/d?' of approximately 0.2 per radian. It will be 
recalled that the special case of a spherical surface with a slit parameter $ = 77 
was solved anal 5 i:ieally in Section IV as a limiting case for a spherical wedge** 
shaped domain. The exact value for the fundamental eigenvalue was analytically 
found to be k= 3/4. The results in Table 8 for $ = 99 tt/ 100, which are included 
in Figure 10, are in excellent agreement with this theoretical result. It is also 
noted from Figure 10, and the results in Table 8 for $ =77 /lOO, that the limiting 
value as $ approaches zero (and the spherical surface with a slit degenerates to 
a sphere with a boimdary point at the north pole) tends rather slowly to the true 
value ^ = 0, with the refinement in the rectangular grids. Further experiments 
conducted with still finer grids in which > 40 indicated that this convergence 
of the fundamental eigenvalue toward zero did indeed proceed, albeit quite 
gradually. 
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Figure 10, The value of the calculated fundamental eigenvalue A. 
as a function of the spherical slit extension parameter $ for 
various rectangular grids (N^ by intervals) 
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APPENDIX A 


ALGORITHM FOR CALCULATING THE FUNDAMENTAL 
EIGENVALUE OF A SPHERICAL TRIANGLE 

1. The followii^ parameters are given initially: 

a. the spherical triangle parameters, 0 , a, b, (7)* 

where 0 < 0 ^ 2 t 7 ^, - oo < a, b < oo; 

b. the rectangular grid parameters, N^ , N^, (13) 

where both N^ and N^ are positive integers greater than imity, 
and, in the general case, must be an even integer; 
e. the parameter, y, defining the grid spacing in the 0 -direction, where (28) 
7 is positive; 

d. the criteria, for iterative convergence, where (41) 

® ^INNER* SuTER ^ 

e. the upper limits, m^a^ * the number of iterations permitted for inner 

and outer convergence, respectively, where both and m„,ax are positive 
integers. 

2. Calculate the grid spacing in the 0 -direction: 

h, = ±. <14a) 

3. Calculate the following vector components as constant parameters for the 
iterative procedure: 

= ih, ^ (13) 

V for i = 1, 2,..., - 1 

M' = a cos 0. + b sin 6, 

1 1 1 J 



for j ^ 0, 1, 2,.. 


-.In^ 


( 28 ) 


♦These equation numbers provide a reference to the discussion in the main body 
of the text. 
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Note: If y = 1, then the above calculations reduce to 


for j = 0, 1, 2, N,, (13; 14b) 

and may be any positive integer greater than unity (not necessarily even). 


N. - cot <h. 
) 


b. r 

j 


-4. ' 

+ 1 - 1 

h2 sin 0. 


2 sin ^ (<h, ^ ct, \ 


c. 

J 


a. - 2b. + c. + c._j 
4. Impose the boundary conditions: 


for i = 0, 1, 2 1 


for j = 1 . 2 ,..., N^- 1 


(24a) 


(24b) 

(24c) 

(23d) 



for j r 0, 1, 2,..., 
for i = 0, 1, 2 


5, Initialize U for the power method: 

if (Mi - N.) s 0, then ^ . = 0; 

i f (Mj - Np < 0, then Uj . = 1 


for i = 1, 2,. 

j - 1, 2, . . 


N, - 1: 
N^- 1. 


(25) 
(27) 

(26) 


6. Initialize the outer iteration index (set m = 1) and initialize the scalar relaxa- 
tion factor (set co = 1) for Gauss-Seidel iteration. 
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7. Calculate: 

*"i.j = for i = 1, 2,..., Kg - 1; 

j - 1. 2,..., N^- 1. (37) 

8. Initialize the inner iteration index (set k= 0). 

9. Increase the value of the inner iteration index by unity, i.e., replace k by (k + 1). 

Also, initialize the convergence parameter by equating it to zero. 


10, Begin the iterative successive overrelaxation process: 
if (M. - N.) ^ 0, then = 0; 

if (M- - Nj) < 0. then 


U 


(k+l) ^ " [b 4- c + b 


.(k) 


(k) 


.(k+l) 






(k) 


for i= 1,2 Ng -1; j=l. 2 - 1. 


(36) 


For each i, j in this range, 

if -uj';>i >€(>'). then set - uj‘;;h 

otherwise, do not change the current value of ef**). This procedure determines 


e<“) = max . (38a) 


11. If k ^2, then calculate the convergence parameter 

rCk) ^ e(k>/g(k-i) (38b) 

In the case k= 1, omit this calculation. 

12. If > Cjpjujgjj and k s then the inner iterations do not converge 

within the allotted number of iterations. If m - 1, disregard such non-convergence 
and proceed to step 13. However, if m > 1, the algorithm does not converge to the 
desired eigenvalue. Return to step 1. 
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If i ^inner then continue the inner iterations by return- 

ing to step 9. 

If e(k) < , then the inner iterations have converged. Proceed to 

step 13. 


13. Calculate the approximate fundamental eigenvalue by the Rayleigh quotient 

‘=1 \j = l ^ y 


(42a) 




\(">) i; P/Q. 


(42b) 


14. If m = 1, proceed to step 15. 

If m > 1 and if - >S”' as in inequality (43), then the outer 

iterations have converged, and is the calculated fimdamental eigenvalue. 
The numerical solution is completed. 

If m > 1 and if ^ Suter inequality (44), then if 

m ^ » the outer iterations do not converge within the allotted number of 

iterations. The algorithm in this case does not converge to the desired eigen- 
value. Return to step 1. 

If m > 1 and if ^ ^outer ♦ if m < m^^^, then continue 

the outer iterations by proceeding to step 15. 


15. Re-initialize U for successive overrelaxation by the power method: 


A » J i r i 


for i = 1, 2,..., - 1; 

j = 1, 2,..., - 1. 


(45) 
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If m = 1 (and assuming that ^ 1), then compute an optimal scalar relaxation 
factor by 

o> = ^ (40) 

1 + /I - rO') 

If m > 1, omit this calculation. 

Increase the value of the outer iteration index by unity, i.e., replace m by 
(m + 1), and resume successive overrelaxation by returning to step 7. 



APPENDIX B 


ALGORITHM FOR CALCULATING THE FUNDAMENTAL 
EIGENVALUE OF A SPHERICAL SURFACE WITH A SUT 

1. The following parameters are given initially: 

a. the slit extension parameter, $, where 0 < § < n', 

b. the rectangular grid parameters, N„, N , 

where both N^ and N^ are positive integers greater than unity, 
and, in the general case, N^ must be an even integer; 

c. the parameter, 7, defining the grid spacing in the 0 -direction, 
where y is positive; 

d. the criteria, , for iterative convergence, 

where 0 < ej^NER . Suter < < and 

e. the upper limits, k , m , on the number of iterations permitted 
for inner and outer convergence, respectively, where both k and 

max 

^raax are positive integers. 

2. Calculate the grid spacing in the (9-direction: 


3. Calculate the following vector components as constant parameters for the 
iterative procedure: 


. 2 iY 77 




i =0, 1, 2 1 n^. 


♦These equation numbers provide a reference to the discussion in the main 
body of the text. 
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Note: If 7 = 1, then the above calculations reduce to 

(13; 54b) 

and N^ may be any positive integer greater than unity (not necessarily even}. 


for j = 0, 1, 2 N^, 


b. = 


(i>. - 0. 


^ T2 — “"7 for i ^ 1, 2,..., N. - 1 

h2 sin 4>. 


2 sinl ( 0 . + 


c. = 

J 


for j =0, 1, 2 N, - 1 


for j = 1 , 2 , - 1 


a. 3 2 b. + c. + c..^ 

4. Impose the boundary conditions: 

Uj = 0 for i = 0, 1, 2 


Uo.i =0 




for j = 0 , 1. 2,.... N^, if 5 


■ 

= 1 | 

0. J 

1 

Ux. 

= 1 


J 


(24a) 

(24b) 

(24c) 

(23d) 

(65) 


(55) 


(56) 


5. Initialize U for the power method: 

Ui.j - 1 for i = 1, 2,..., - 1: j = l, 2.....N, 


for j = 1 , 2 ,.... N^, if >$. 


6. Initialize the outer iteration index (set m = 1) and initialize the scalar relax- 
ation factor (set w = 1) for Gauss-Seidel iteration. 


7. Calculate: 

i = 0, 1, 2,.... N^; j = l, 2,...,N, 




(69) 
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8. Initialize the inner iteration index (set k= 0). 

9. Increase the value of the inner iteration index by unity, i.e., replace k by 

(k + 1). Also, initialize the convergence parameter by equating it to 

zero. 

10. Begin the iterative successive overrelaxation process: 




l(l') 


(k) 


,(k-+l) 


+ c 


+ +<■!-<») u;';'. 


(k) 


for i = 1, 2 - 1; j = 1, 2,..., - 1. 

FoUowing the above computation for i * N - 1; j, compute: 


(36) 




f o 


r j - 1, 2, . . . , - 1, if > $; 


(67) 

(57) 


f°r J = 1, 2 N^-1, if0, >$. 

Now proceed with the computations (36), interspersing the computations (67) and 
(57) as necessary. 

Upon completion of the computations (36), compute: 

F. 


U 


(k+l) _ (U) 


1 . N 


4> N 


'0 


T uf*‘> . ±a*. 
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for i= 0, 1, 2,..., Ng. 

For each i= 0, 1, 2, . . . , N^; j = 0, 1, 2 N^, if | 

then set ^ “ I U| ^ ^ - Uj*^ • {; otherwise, do not change the current value of 
^ . This procedure determines 

eC') = max I - 


0 < i < N 


'0 


.0 < j < N 


0 


(70) 
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11. If k i 2, then calculate the convergence parameter 

r(k) ^ g(k)/^(k-i)_ (3gbj 


In the case k= 1, omit this calculation. 

12. If € > and k then the inner iterations do not converge 

within the allotted number of iterations. If m = 1, disregard such non- 
convergence and proceed to step 13. However, if m > 1, the algorithm does not 
converge to the desired eigenvalue. Return to step 1. 

If and k < k^^^, then continue the inner iterations by return- 

ing to step 9. 

If €(‘'> < then the inner iterations have converged. Proceed to step 13. 

13. Calculate the approximate fundamental eigenvalue by the Rayleigh quotient: 



Q = 



F. .U. 

i . j 1 . j 


+ F. .. U, 




= P/Q. 


(71b) 


14. If m = 1, proceed to step 15. 

If m > 1 and if IxC") - as in inequality (43), then the outer 

iterations have converged, and is the calculated fundamental eigenvalue. 

The numerical solution is completed. 

If m > 1 and if ^ ^oute» as in inequality (44), then if 

™ - ^niax » outer iterations do not converge within the allotted number of 
iterations. The algorithm in this case does not converge to the desired eigenvalue. 
Return to step 1. 

If m > 1 and if |\<"’) - XC"--!)! > , but if m < m„^^. then continue the 

outer iterations by proceeding to step 15. 
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If m >1, omit this calculation. 

Increase the value of the outer iteration index by unity, i.e., replace m by 
(m + 1), and resume successive overrelaxation by returning to step 7, 



